# R script for paper Manger - Pickup
# The Coevolution of Trade Agreement Networks and Democracy
# Journal of Conflict Resolution

#setwd("pathname /Replication Files")

# read data
net83 <- as.matrix(read.table("CoevPTA1983.dat"))
net85 <- as.matrix(read.table("CoevPTA1985.dat"))
net86 <- as.matrix(read.table("CoevPTA1986.dat"))
net88 <- as.matrix(read.table("CoevPTA1988.dat"))
net91 <- as.matrix(read.table("CoevPTA1991.dat"))
net92 <- as.matrix(read.table("CoevPTA1992.dat"))
net93 <- as.matrix(read.table("CoevPTA1993.dat"))

democ  <- as.matrix(read.table("udsInt83-93_1.dat"))
democ1<-cbind(democ[1:145,1],democ[1:145,3],democ[1:145,4],democ[1:145,6],democ[1:145,9],democ[1:145,10],democ[1:145,11])
invGDP <- as.matrix(read.table("invGDP83-93.dat"))
invGDP1<-cbind(invGDP[1:145,1],invGDP[1:145,3],invGDP[1:145,4],invGDP[1:145,6],invGDP[1:145,9],invGDP[1:145,10],invGDP[1:145,11])
GDPcap <- as.matrix(read.table("GDPcap83-93.dat"))
GDPcap1<-cbind(GDPcap[1:145,1],GDPcap[1:145,3],GDPcap[1:145,4],GDPcap[1:145,6],GDPcap[1:145,9],GDPcap[1:145,10],GDPcap[1:145,11])

dist <- as.matrix(read.table("Coevdist.dat"))

Ally83 <- as.matrix(read.table("CovAlly1983.dat"))
Ally85 <- as.matrix(read.table("CovAlly1985.dat"))
Ally86 <- as.matrix(read.table("CovAlly1986.dat"))
Ally88 <- as.matrix(read.table("CovAlly1988.dat"))
Ally91 <- as.matrix(read.table("CovAlly1991.dat"))
Ally92 <- as.matrix(read.table("CovAlly1992.dat"))


Trade83 <- as.matrix(read.table("Cvtrade1983.dat"))
Trade85 <- as.matrix(read.table("Cvtrade1985.dat"))
Trade86 <- as.matrix(read.table("Cvtrade1986.dat"))
Trade88 <- as.matrix(read.table("Cvtrade1988.dat"))
Trade91 <- as.matrix(read.table("Cvtrade1991.dat"))
Trade92 <- as.matrix(read.table("Cvtrade1992.dat"))

# Prepare data for RSiena.
library(RSiena)

# create dependent variable
PTAs <- sienaNet(array(c(net83, net85, net86,
               net88, net91, net92, net93),
                        dim = c(145, 145, 7) ))

#create dependent behavior variable 
Democracy <- sienaNet(democ1, type="behavior")

# create monadic time varying covariates
GDP_inv   <- varCovar(invGDP1)

# create constant dyadic covariate
Distance  <- coDyadCovar(dist)

# create dyadic time varying covariates
Alliance  <- varDyadCovar(array(c(Ally83, Ally85, Ally86,
               Ally88, Ally91, Ally92),
                          dim = c(145, 145, 6) ))

Trade   <- varDyadCovar(array( c(Trade83, Trade85, Trade86,
               Trade88, Trade91, Trade92),
                        dim = c(145, 145, 6) ))


# define data                         
PTAbehdata <- sienaDataCreate(PTAs, Democracy, GDP_inv, Distance,
                           Alliance, Trade)
# create effects structure
PTAeff <- getEffects(PTAbehdata )

# give global report about data
print01Report(PTAbehdata, PTAeff, modelname = 'PTA_coev8393_uds_int' )
   
# Define the kind of non-directed model 
model33 <- sienaModelCreate(projname = "PTA_coev8393_uds_int", modelType = 3, n3=3000)

# Specify the model
PTAeff <- setEffect(PTAeff, transTriads, include=FALSE)
PTAeff <- includeEffects(PTAeff, nbrDist2)
PTAeff <- includeEffects(PTAeff, X, interaction1="Distance")
PTAeff <- includeEffects(PTAeff, X, interaction1="Alliance")
PTAeff <- includeEffects(PTAeff, altX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, altX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, X, interaction1="Trade")
PTAeff <- includeInteraction(PTAeff, X, altX, interaction1=c("Trade", "GDP_inv"))
PTAeff <- includeEffects(PTAeff, name = "Democracy", avAlt, outdeg, interaction1 = "PTAs")
PTAeff <- setEffect(PTAeff, name = "Democracy", quad, include=TRUE)

# Check the current model specification
PTAeff

# estimate
ans2 <- siena07(model33, data = PTAbehdata, effects = PTAeff)
ans2

# continue estimate
ans2 <- siena07(model33, data = PTAbehdata, effects = PTAeff, prevAns=ans2)
ans2

# continue estimate
ans2 <- siena07(model33, data = PTAbehdata, effects = PTAeff, prevAns=ans2)
ans2

# continue estimate
ans2 <- siena07(model33, data = PTAbehdata, effects = PTAeff, prevAns=ans2)
ans2